function [grad,hess1,hess2] = Full_HessScore(x,f)
% This file reads X, points of treamtent, and F, the vectors of scores, to
% create two hessian variants and the gradient. X should be (2*k x k) in
% dimension, consisting of symmetric points around the minimizing
% parameters. F should be (obs x 2*k), giving the vector of points for each
% run of X.
%
% Hess1 is B in langauge of Train (the average outer product of scores),
% while Hess2 is W (the covariance of scores around their mean). When at
% the true minimum, the gradient is zero and W=B.

obs = size(f,1);
k   = size(f,2)/2;

%% Step 1: Calculate scores (s)

s   = nan(obs,k); 

for j=1:k
    s(:,j) = (f(:,2*j) - f(:,2*j-1)) ./ (2*(x(2*j-1,j)-x(2*j,j)));
end

%% Step 2: Calculate gradient

grad = sum(s,1)/obs;

%% Step 3: Calculate Hess1

hessop = zeros(k,k);

for i=1:obs
    hessop = hessop + s(i,:)'*s(i,:);
end

hess1 = hessop/obs;

%% Step 4: Calculate Hess2

hesscov = zeros(k,k);

for i=1:obs
    hesscov = hesscov + (s(i,:)-grad)'*(s(i,:)-grad);
end

hess2 = hesscov/obs;

end